In kongdd/PhenoAsync: Extract Remote Sensing Vegetation Phenology

1. select available sites first

• (1) Site-year with 300 valid daily observations were left.
• (2) Days with GPP_NT >= 1 should be more than 100.

Considering to remove those sites: * Remove bad sites US-Tw3 : only 1 year data, 5 growing season in it. CH-Fru : has checked NDVI, EVI and GPP_mod, only one growing season. AU-Dyu : too many missing values, suggest to remove this site AU-Cum : only 1 year data, not seasonality AU-GWW : only 1 year data, not seasonality AU-Wac : only 1y, and VI is inverse of GPP

GF-Guy : EBF, not seasonality AU-Cum : EBF, not seasonality AU-Cpr : negative corr with VI

source('inst/shiny/check_season/global.R')
source("test/stable/load_pkgs.R")
# source("test/phenology_async/R/s1_materials/main_phenofit.R")
load("data/phenoflux_115_gs.rda")

st[, `:=`(IGBPname = IGBP, lon = long)]

check grow season dividing

1. Check check_input function
2. Check wWHIT performance

df <- lst_sm$MOD13A1
i  <- 1

# source('inst/shiny/check_season/global.R')
sites_part <- st$site # [IGBP == "CRO"]
sites <- sites_part
IGBP_forest <- c("DBF", "EBF", "ENF", "MF", "CSH", "GRA")
IGBP_whit   <- c("CRO")

nptperyear <- 365
file <- sprintf("fluxsite GPP_DT Growing season dividing v0.2.0.pdf")
CairoPDF(file, width = 10, height = 2*6)
par(mfrow = c(6, 1), mar = c(2, 3, 2, 1), mgp = c(1.5, 0.6, 0))

# sites <- unique(d_obs$site)
res   <- list()
stats <- list()
for (i in seq_along(sites)){
    runningId(i)
    sitename <- sites[i]
    tryCatch({
        # check_season(sitename, df_raw, stations)
        d     <- df[site == sitename & scale == "0m", .(t, date, y = EVI, w, SummaryQA)] #%T>% plotdata(365)
        sp    <- st[site == sitename, ]

        south <- sp$lat < 0 
        d_new <- add_HeadTail(d, south = south)

        nptperyear <- 23
        INPUT <- with(dnew, check_input(t, y, w, nptperyear, maxgap = ceiling(nptperyear/12*1.5)))
        INPUT$south <- sp$lat < 0
        plotdata(INPUT)

        # parameters for season_3y
        threshold_max = 0.1
        nf = 1

        FUN_fit       <- "wWHIT"
        threshold_max <- ifelse(cv_coef(d$y)[3] >= 1, 0.1, 0.2) # experience param
        # FUN_fit <- ifelse(sp$IGBP %in% IGBP_forest, "wHANTS", "wWHIT")
        wFUN <- wTSM# "wBisquare"
        maxExtendMonth <- ifelse(sp$IGBP == "EBF", 2, 2)

        # wFUN <- "wBisquare", "wTSM", threshold_max = 0.1, IGBP = CSH
        # INPUT <- get_input(df, st, sitename)
        brks2  <- season_3y(INPUT, south = INPUT$south, rFUN = get(FUN_fit),
                       wFUN = wFUN,
                     IsPlot = IsPlot,
                     lambda = 10,
                     iters = 2,
                     minpeakdistance = nptperyear/6,
                     MaxPeaksPerYear = 3,
                     MaxTroughsPerYear = 4,
                     ypeak_min = 0.08, 
                     IsOnlyPlotbad = FALSE
        ## Get curve fitting
    )

    if (IsPlot){
        abline(h = 1, col = "red")
        title(INPUT$titlestr)
    }

        res[[i]] <- brks 
    })
}

dev.off()
names(res) <- sites
brks_lst <- res
# save(brks_lst, file = "data_test/phenoflux_115_gs.rda")

kongdd/PhenoAsync documentation built on April 21, 2024, 2:36 a.m.